function [Tinert2crtbp,Tcrtbp2inert,S_eci] = CRTBP2Inertial(JD_interest,S_crtbp,JD, R2, V2,body)
%% Constants
G           = 6.674e-11 * ((1/1000)^3);
mass_Mars   = 6.39e23;
mass_Deimos = 1.8e15;
mass_Phobos = 10.8e15;
GM_Mars     = G*mass_Mars;
GM_Deimos   = G*mass_Deimos;
GM_Phobos   = G*mass_Phobos;
GM_sun      = 1.327124400*10^(11);
%[km^3/s^2] Gravitational Parameters

u_sm = GM_Mars / (GM_Mars+GM_sun);
u_md = GM_Deimos / (GM_Deimos+GM_Mars);
u_mp = GM_Phobos / (GM_Phobos+GM_Mars);
%[] Mu constant
%    sm --> Sun-Mars system
%    em --> Earth-Moon system%%

DU_sm = 235.34e6;
DU_mp = (9236.53+9517.58)/2;
DU_md = (23455.5+23470.9)/2;
%[km] Distance from Sun to Mars

TU_sm = sqrt(DU_sm^3/GM_sun);
TU_mp = sqrt(DU_mp^3/GM_Mars);
TU_md = sqrt(DU_md^3/GM_Mars);
%[s] Time unit
%[] Parameter for altering the halo orbit's state

%%

if body == 'Phobos'
    GM_1= GM_Mars;          %Earth
    GM_2 = GM_Phobos;       %Phobos
    TU = TU_mp;
    LU = DU_mp;
    VU = LU/TU;
    %[] Canonical unit conversions
    Wcm = sqrt((GM_1+GM_2)/LU^3)*[0;0;1];
    %[] angular velocity in Cr3bp frame
elseif body == 'Deimos'
    GM_1= GM_Mars;          %Mars
    GM_2 = GM_Deimos;       %Deimos
    TU = TU_md;
    LU = DU_md;
    VU = LU/TU;
    %[] Canonical unit conversions
    Wcm = sqrt((GM_1+GM_2)/LU^3)*[0;0;1];
    %[] angular velocity in Cr3bp frame
end
%%



S_crtbp = [S_crtbp(1:3)*LU;S_crtbp(4:6)*VU];
R_crtbp = S_crtbp(1:3);
V_crtbp = S_crtbp(4:6);

%[] Non-inertial position and velocity of the spacecraft at JD_interest.
%% 1: Using julian date, determine moon's state at interested julian date.

JD_difference_vector = JD-JD_interest;
%[] Determine the difference time vector

ind = find(JD_difference_vector > 0 , 1 , 'first') - 1;
%[] Find the index of most close previous state

R2_o = R2(:,ind);
V2_o = V2(:,ind);
S2_o = [R2_o;V2_o];
%[km,km/s] Position and velocity of the moon WRT Earth.

UTC_JD_interest = CalendarDate(JD_interest);
UTC_JD_closest = CalendarDate(JD(ind));
%[] Get UTC time of JD time of interest, which is the JD_interest and
%JD(ind), which is the closest previous information.

time_difference_seconds = etime(UTC_JD_interest,UTC_JD_closest);

%[] Determine time difference in seconds.
if time_difference_seconds == 0
    R2_interest = R2_o;
    V2_interest = V2_o;
    S2_interest = [R2_interest;V2_interest];
else
    [~,S] = ode45(@(t,S)R2bpCartesianModel(t,S,GM_1),[0,time_difference_seconds],S2_o);
    S = S';
    %[] Integrate for difference amount to determine exact state of the moon
    %wrt Earth in ECI frame.
    
    R2_interest = S(1:3,end);
    V2_interest = S(4:6,end);
    S2_interest = [R2_interest;V2_interest];
    %[] Position and velocity of the moon at JD_interest.
end
%% 2: Determine states of CM

Scm1 = (GM_2*S2_interest)/(GM_2+GM_1);
Rcm1 = Scm1(1:3);
% Vcme = Scme(4:6);
%[] Determine the state of the center of mass at JD of interest
%There are no additional terms because the Earth's position wrt Earth is
%zero vector and the gravitational parameter of the satellite is assumed to
%be zero.

%% 3: Make all states WRT CM

S1cm = -Scm1;
% Rearthcm = Searthcm(1:3);
% Vearthcm = Searthcm(4:6);
%[] State of the Earth wrt CM

S2cm = S2_interest-Scm1;
% Rmooncm = Smooncm(1:3);
% Vmooncm = Smooncm(4:6);
%[] State of the moon wrt CM

% Ssatcm = S_eci-Scme;
% Rsatcm = Ssatcm(1:3);
% Vsatcm = Ssatcm(4:6);
%[] State of the state of interest [S_eci] wrt Center of mass.

%% 4: Determine Teci2crtbp and Tcrtbp2eci

Hm          = cross(R2_interest,V2_interest);
%[1/s]Current specific angular momentum of the Moon WRT the CM in ECI coordinates.

Rhat        = R2_interest / norm(R2_interest);
%[]Current CRTBP radial direction in ECI coordinates.

Nhat        = Hm / norm(Hm);
%[]Current CRTBP normal direction in ECI coordinates.

That        = cross(Nhat,Rhat);
%[]Current CRTBP tangential direction in ECI coordinates.

Tcrtbp2inert  = [Rhat, That, Nhat];
%[]Matrix that transforms vectors from CRTBP to ECI coordinates.

Tinert2crtbp  = transpose(Tcrtbp2inert);
%[] Matrix that transforms vectors from ECI to CRTBP coordinates.


%% 5: Determine S_crtbp

Rsatcm_eci  = Tcrtbp2inert * R_crtbp;
Vsatcm_eci  = Tcrtbp2inert*(V_crtbp+cross(Wcm,Tinert2crtbp*Rcm1)+cross(Wcm,R_crtbp));
Ssatcm_eci  = [Rsatcm_eci;Vsatcm_eci];
S_eci       = Ssatcm_eci + Scm1;

end